function x = suprarelaxare1(A, b, x0, w, eps, maxiter)
%	Aceasta functie rezolva sistemul A*x = b; Pentru a se rezolva
%acest sistem mai usor, se considera A = N - P; (N - P)*x = b;
%N*x = P * x + b; fie N1 inversa lui N; x = N1*P*x + N1*b. Constructia 
%matricilor N si P pleaca de la matricile D, L, U. D este o matrice
%formata din diagonala principala a matricei A, si are 0 in afara 
%diagonalei principale. L este o matrice inferior triunghiulara, 
%care are sub diagonala principala elementele lui A de sub diagonala 
%principala( luate cu semn schimbat) si 0 in rest. U este o matrice
%superior triunghiulara; are deasupra diagonalei principale elementele 
%lui A de deasupra diagonalei principale( luate cu semn schimbat) si 
% 0 in rest. In aceasta metoda, N si P depind de parametrul w, numit 
%parametru de relaxare. Astfel, N = 1/w * D - L si P = (1/w - 1) * D + U .
%Pornind de la relatia: x = N1*P*x + N1*b, solutia se calculeaza recursiv.
%Astfel, x(la pasul curent) = N1*P*x( la pasul anterior) + N1*b. De aceea,
%se pleaca de la un x0( o aproximare initiala) si se calculeaza o alta 
%aproximare, pana cand || A * x - b|| < eps
%sau pana cand numarul de iteratii depaseste numarul maxiter.
%
%	Date de intrare:
%		-	A  ->  marricea cu care se inmulteste necunoscuta in sistemul  A*x = b;
%		-	b  ->  coloana termenilor liberi;
%		-	x0 ->  o aproximare initiala a solutiei;
%		-	w  ->  parametru de suprarelaxare;
%		-	eps->  toleranta cu care se calculeaza solutia sistemului;
%		-maxiter-> numarul maxim de iteratii.
%
%	Date de iesire:
%		-	x  -> solutia sistemului.

n = length(b); %aflam dimensiunea matricei A

%initializam matricile D, L, U cu O
D = zeros(n, n);
L = zeros(n, n);
U = zeros(n, n); 

%calculam efectiv cele 3 matrici

%matricea D este formata din diagonala matricei A si zero in rest
for i = 1:n
	D(i, i) = A(i, i);
endfor

%matricea D este formata din elementele de deasupra diagonalei principale ale matricei A(luate cu semnul -)  si o in rest
for i = 1:n
	for j = i+1:n
		U(i, j) = -A(i, j);
	endfor
endfor

%matricea L este formata din elementele de sub diagonala principala a matricei A, luate cu semnul - si 0 in rest
for i = 1:n
	for j = 1:i - 1
		L(i, j) = -A(i, j);
	endfor
endfor

%acuma calculam matricile N(w) si P(w), dupa urmatoarele formule:
%N(w) = 1/w * D - L;
%P(w) = (1/w - 1) * D + U
for i = 1:n
	for j = 1:n
		N(i, j) = (1/w)*D(i, j) - L(i, j);
		P(i, j) = (1/w - 1)*D(i, j) + U(i, j);
	endfor
endfor

%calculam separat N^(-1) si retinem rezultatul in N1
N1 = inv(N);

%calculam matricile G si c
G = N1 * P;
c = N1 * b;

%incepem sa calculam pe x; pentru aceasta, retinem in "iter" numarul de iteratii de la
%un moment-dat; iter <= maxiter

iter = 0;
x = G * x0 + c;%se calculeaza un x in afara ciclului 

while( iter <= maxiter  )
	iter = iter + 1;
	x = G * x0 + c;
	if( norm( x - x0) < eps) 
		break;
	endif
	x0 = x;%in x0 se retine acel x calculat la pasul anterior\
	
endwhile

endfunction
